In this file, the goal is to quantify the historical cloud cover in the planned flight boxes and to assess the implications for BioSCape operations.
# Load required packages
# Load packages
library(rgee)
library(targets)
library(sf)
library(tidyverse)
library(lubridate)
library(leaflet)
library(ggplot2)
library(leafem)
#Load required data
# get domain
domain <- st_read("data/output/domain.gpkg",
quiet = TRUE)
domain_sf <- domain
# get flight boxes
boxes <- st_read("data/flight_planning/20221026_flightboxes.gpkg",
quiet = TRUE)
boxes$id <- 1:20 # need a unique ID to make things easier
boxes_sf <- boxes
# Download table from drive (to see the code underlying this or to update the data, see the file "R/mock_flights_earth_engine.R")
googledrive::drive_download(file = "EMMA/cloud_stats.csv",
path = "data/test_cloud_stats.csv",
overwrite = TRUE) #note: if dling more than one need to add a prefix or something
cloud_table <- read.csv("data/test_cloud_stats.csv")
# Parse date
cloud_table %>%
mutate(year = year(date),
month = month(date),
day = day(date),
day_of_year = yday(date)) -> cloud_table
ggplot()+
geom_sf(data = domain)+
geom_sf(data = boxes_sf,fill="light green")+
geom_sf_text(data = boxes_sf,aes(label = id),size=10)+
labs(x=NULL,y=NULL)
Figure 1. The full domain (light grey) and planned flight boxes (green). Numbers shown are box IDs (used below).
Focal study sites have been selected on the basis of the following criteria:
#ID by day of year
cloud_table %>%
group_by(id,day_of_year)%>%
ggplot() +
geom_tile(mapping = aes(x = day_of_year,
y = id,
fill = mean))+
scale_fill_gradient(low = "sky blue",
high = "white")+
facet_wrap(~year)+
labs(fill="Mean\nCloud\nCover",
y="ID",
x = "Day of Year")
Figure 2. Rows represent different flight boxes (see Fig 1.), columns different days.
# mean monthly cloud cover
cloud_table %>%
group_by(id, month)%>%
filter(month != 9)%>%
summarize(mean_cc = mean(na.omit(mean)))%>%
inner_join(x = boxes_sf)%>%
ggplot(mapping = aes(fill = mean_cc))+
geom_sf()+
geom_sf(data = domain_sf,
inherit.aes = FALSE,fill=NA)+
scale_fill_gradient(low = "sky blue",high = "white")+
facet_wrap(~month,ncol = 1)+
geom_sf_text(aes(label = round(mean_cc,digits = 2)))+
labs(fill="Mean\nCloud\nCover")+
xlab(NULL)+
ylab(NULL)
Figure 3. The full domain (light grey) and planned flight boxes, colored by mean monthly cloud cover.Numbers on flight boxes refer to mean cloud cover.
Here, we define a “clear” day as one with a mean cloud cover of less than 10 percent.
#prop clear days
cloud_table %>%
na.omit()%>%
filter(month != 9)%>%
mutate(binary_clear = dplyr::if_else(mean <= .1,true = 1,false = 0)) %>%
group_by(id)%>%
summarize(prop_clear = sum(binary_clear)/n(),
clear_days = sum(binary_clear),
total_days = n())%>%
inner_join(x = boxes_sf)%>%
ggplot(mapping = aes(fill = prop_clear))+
geom_sf()+
geom_sf(data = domain,inherit.aes = FALSE,fill=NA)+
scale_fill_gradient(low = "white",high = "sky blue",limits=c(0,1))+
geom_sf_text(aes(label = round(prop_clear,digits = 2)),size=10)+
labs(fill = "Prop.\nClear",
x=NULL,
y=NULL)
Figure 4. Proportion of clear days (mean cloud cover < 10% ) for each flight box. Values in boxes are the proportion of clear days for that box
The map below is interactive and allows zooming, moving, and toggling layers on/off.
#Pull in other bioscape layers
cloud_table %>%
na.omit()%>%
filter(month != 9)%>%
mutate(binary_clear = dplyr::if_else(mean <= .1,true = 1,false = 0)) %>%
group_by(id)%>%
summarize(prop_clear = sum(binary_clear)/n(),
clear_days = sum(binary_clear),
total_days = n())%>%
mutate(prop_clear = round(prop_clear,2))%>%
inner_join(x = boxes_sf)%>%
st_transform(crs = st_crs(4326)) -> flights
team_requests <- st_read("data/manual_downloads/BIOSCAPE_proposed/20221014_team_polygons.gpkg",quiet = TRUE) %>%
st_transform(crs = st_crs(4326))
domain_sf %>%
st_transform(crs = st_crs(4326)) -> domain_wgs84
#Make a palette
pal <- colorNumeric(palette = colorRamp(c("white", "blue"), interpolate = "spline"),
domain = unique(flights$prop_clear))
#Make labels
labels <- sprintf(as.character(flights$prop_clear)) %>%
lapply(htmltools::HTML)
leaflet(data = flights) %>%
addProviderTiles("Esri.NatGeoWorldMap", group = "NatGeo") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "World Imagery") %>%
addMapPane("flights", zIndex = 420) %>%
addMapPane("requests", zIndex = 410)%>%
addPolygons(stroke = TRUE,
group = "Flights",
color = ~pal(prop_clear),
opacity = 1,
label = labels,
options = pathOptions(pane = "flights"))%>%
addPolygons(data = domain_wgs84,
stroke = TRUE,
color = "grey",
fill = FALSE,
weight = 3) %>%
addPolygons(data = team_requests%>%
st_zm(drop = T, what = "ZM"),
stroke = TRUE,
color = "brown",
group = "Requests",
options = pathOptions(pane = "requests"))%>%
addMouseCoordinates() %>%
#addImageQuery(sampling_options_wgs84, type="mousemove", layerId = "park_name") %>%
leaflet::addLegend(position = "bottomright",
pal = pal,
values = unique(flights$prop_clear),
opacity = 1,
title = "Cluster") %>%
addLayersControl(
baseGroups = c("World Imagery","NatGeo"),
overlayGroups = c("Flights","Requests"),
options = layersControlOptions(collapsed = FALSE),
position = "topright")